function Q = fourier(x, fs, f0)
% 傅里叶系数 Q
%
% 输入:
%   x  : 神经元膜电位的时间序列 (行或列向量)
%   fs : 采样频率 (Hz)，若仿真步长为 h，则 fs = 1/h
%   f0 : 目标频率 (Hz)，即想要计算 Q 的频率
%
% 输出:
%   Q_val : 目标频率 f0 对应的傅里叶系数 Q

    x = x(:);                % 转为列向量，防止维度不一致
    x_dc = x - mean(x);      % 去除直流(平均值)

    % 计算快速傅里叶变换 (FFT)
    N = length(x_dc);
    X = fft(x_dc);
    X = X(1:floor(N/2)+1);      % 仅保留正频率部分（含直流）

    % 计算频率分辨率并确定目标频率索引
    freqResolution = fs / N;
    kf = round(f0 / freqResolution);

    % 检查目标频率 f0 是否在有效范围 [0, fs/2]
    if f0 < 0 || f0 > fs/2
        error('目标频率 f0 必须在 [0, %.2f] Hz 范围内。', fs/2);
    end

    % 取幅值并做归一化
    % abs(X(kf+1)) 是目标频率下的复数幅值模
    Q = abs(X(kf+1)) / N;
end
